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We present the first calculation of gravitational wave emission produced in the gravitational collapse of uni¬ 
formly rotating neutron stars to black holes in fully three-dimensional simulations. The initial stellar models are 
relativistic polytropes which are dynamically unstable and with angular velocities ranging from slow rotation to 
the mass-shedding limit. An essential aspect of these simulations is the use of progressive mesh-refinement tech¬ 
niques which allow to move the outer boundaries of the computational domain to regions where gravitational 
radiation attains its asymptotic form. The waveforms have been extracted using a gauge-invariant approach in 
which the numerical spacetime is matched with the non-spherical perturbations of a Schwarzschild spacetime. 
Overall, the results indicate that the waveforms have features related to the properties of the initial stellar mod¬ 
els (in terms of their w-mode oscillations) and of the newly produced rotating black holes (in terms of their 
quasi-normal modes). While our waveforms are in good qualitative agreement with those computed by Stark 
and Piran in two-dimensional simulations |.1J, our amplitudes are about one order of magnitude smaller and 
this difference is mostly likely due to our less severe pressure reduction. For a neutron star rotating uniformly 
near mass-shedding and collapsing at 10 kpc, the signal-to-noise ratio computed uniquely from the burst is 
S/N ~ 0.25, but this grows to be S/N < 4 in the case of LIGO II. 

PACS numbers: 04.25.Dm, 04.30.Db, 04.70.Bw, 95.30.Lz, 97.60.Jd 


The study of the gravitational collapse of rotating stars to 
black holes is a cornerstone of any theory of gravity and a 
long standing problem in general relativity. Over the years, 
numerous approaches have been applied and several different 
techniques developed to tackle this problem which is not just 
academic. Indeed, important issues in relativistic astrophysics 
awaiting clarification, such as the mechanism responsible for 
7-ray bursts, may be unveiled with a more detailed under¬ 
standing of the physics of gravitational collapse in rotating 
and magnetized stars. Furthermore, the study of gravitational 
collapse will provide the waveforms and the energetics of one 
of the most important sources of gravitational radiation. 

In the absence of symmetries, no analytic solutions are 
known that describe the gravitational collapse to a black hole 
and it is therefore through numerical relativity simulations 
that one can hope to improve our knowledge of this pro¬ 
cess under realistic conditions. However, this is not an easy 
task and the modelling of black hole spacetimes with col¬ 
lapsing matter-sources in multidimensions is one of the most 
formidable efforts of numerical relativity. This is due both 
to the inherent difficulties and complexities of the system of 
equations which is to be solved ( i.e. the Einstein equations 
coupled to the relativistic hydrodynamics equations) and to 
the vast computational resources needed in three-dimensional 
(3D) evolutions. 

In addition to the technical difficulties due to the accurate 
treatment of the hydrodynamics involved in the collapse, the 
precise calculation of the gravitational radiation emitted in the 
process is particularly challenging as the energy released in 
gravitational waves is much smaller than the total rest-mass 
energy of the system. Indications of the difficulties inher¬ 
ent to the problem of calculating the gravitational-wave emis¬ 


sion in rotating gravitational collapse have emerged in the first 
and only work on this, which dates back to almost 20 years 
ago d. In 1985, in fact, in a landmark work in numerical 
relativity. Stark and Piran used their axisymmetric general- 
relativistic code to evolve rotating configurations and to com¬ 
pute the gravitational radiation produced by their collapse to 
black holes. The results referred to initial configurations con¬ 
sisting of polytropic stars which underwent collapse after the 
pressure was reduced by a factor ranging from 60% up to 99% 
for the rapidly rotating models. The initial data effectively 
consisted of spherically-symmetric solutions with a uniform 
rotation simply “added” on. Not being stationary solutions of 
the Einstein equations, these stars could reach dimensionless 
spins up to a = J/M 2 = 0.94, with J and M the angular 
momentum and mass of the star, respectively. 

Overall, their investigation revealed that while the nature 
of the collapse depended on the parameter a, the form of the 
waves remained roughly the same over the entire range of the 
values of a, with the amplitude increasing with a. Particularly 
important was the evidence that, despite the complex dynam¬ 
ics of the matter during the collapse, the gravitational-wave 
emission could essentially be related to the oscillations of a 
perturbed black hole spacetime. 

In recent years many studies have extended to three spa¬ 
tial dimensions the investigation of gravitational collapse to 
black holes OS01- Despite the improvements in the evolu¬ 
tion of the hydrodynamics and of the spacetime achieved by 
these simulations, none of them has addressed the problem 
of the gravitational-wave emission. The reason for this is to 
be found in the small amplitude of the signal, often below the 
truncation error of the 3D simulations, but most importantly in 
the fact that all of the above calculations made use of Cartesian 
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grids with uniform spacing. With the computational resources 
currently available, this choice places the outer boundaries too 
close to the source to detect gravitational radiation. 

Using a recently developed code for the solution of the Ein¬ 
stein equations in non-vacuum spacetimes, the Whisky code, 
we have investigated the collapse of rapidly rotating relativis¬ 
tic stars to Kerr black holes y|]. An important aspect of these 
simulations is the detailed study of the geometrical and dy¬ 
namical properties of both the apparent and event horizons, 
allowing for the determination of both the mass and angular 
momentum of the black hole with an accuracy unprecedented 
for a 3D simulation. In turn, these measures have set upper 
limits on the energy and angular momentum lost during the 
collapse in the form of gravitational radiation, the first such 
estimates coming from 3D calculations. 

However, as in previous works (e.g. refs. E 00 X the sim¬ 
ulations reported in 0 made use of numerical grids with uni¬ 
form spacing and thus with outer boundaries very close to the 
initial position of the stellar surface. With these restrictions, 
the gravitational radiation extracted does not provide interest¬ 
ing information besides the obvious change in the quadrupole 
moment of the background spacetime. As we will show 
below, the use of progressive mesh-refinement (PMR) tech¬ 
niques has removed these restrictions, enabling us to place the 
outer boundaries of the computational domain at very large 
distances from the collapsing star. This conceptually simple 
but technologically challenging improvement has two impor¬ 
tant physical consequences. Firstly, it reduces the influence of 
inaccurate boundary conditions at the outer boundaries of the 
domain whilst retaining the required accuracy in the region 
where the black hole forms. Secondly, it allows the wave- 
zone to be included in the computational domain and thus to 
extract the gravitational waves produced in the collapse. 

In practice, we have adopted a Berger-Oliger prescription 
for the refinement of meshes on different levels jg] and used 
the numerical infrastructure described in 0 . In addition to 
this, we have also implemented a simplified form of adaptivity 
in which new refined levels are added at predefined positions 
during the evolution. More specifically, given an initial stel¬ 
lar model of mass M and equatorial coordinate radius R, our 
initial grid consists of four levels of refinement, with the inner¬ 
most one covering the star with a typical resolution of Ax ~ 
0.17 M and with the outermost having a typical resolution of 
Ax ~ 1.38 M and extending up to ~ 82.5 M ~ 20.9 R. As 
the collapse proceeds and the star occupies smaller portions of 
the computational domain, three more refined levels are added 
one by one, nested in the four original ones. By the time the 
simulation is terminated at ~ 81.5 M, the finest typical spatial 
resolution is Ax ~ 0.02 M. A detailed discussion of the grid 
and of its evolution will be given in 0 . 

The initial data for our simulations is the same described in 
iBll and basically consists of axisymmetric rotating relativistic 
stars, calculated as equilibrium solutions of the Einstein equa¬ 
tions in a compactified domain and in polar coordinates. For a 
direct comparison with the results in Jlj] and because no shock 
is observed during the collapse, the stars are modelled with a 


polytropic equation of state (EOS) p = Kp r , with T = 2 and 
with the polytropic constant which is initially /i ID = 100 . 
Once secularly unstable solutions are found along sequences 
of fixed angular momentum or fixed rest-mass, we consider 
initial models that have the same axis ratios but slightly larger 
central energy densities and are dynamically unstable. Here¬ 
after we will restrict the discussion to the collapse of the most 
rapidly rotating dynamically unstable model, namely model 
D4 in |0, which represents a star of mass M = 1.861 M 0 , 
circumferential equatorial radius R e = 14.25 km and rotating 
close to the mass-shedding limit with a = 0.54 (cf Table I 
of 0 ). The discussion of the emission from stellar models 
rotating at smaller velocities will be presented in 0 . 

Although model D4 is dynamically unstable, it is very close 
to a stationary solution and we trigger its collapse by reduc¬ 
ing the pressure by < 2 %, with a change in the value of the 
polytropic constant. After the perturbation is introduced, the 
constraint equations are again solved to enforce the constraint 
violation to be at the truncation-error level. As a contrast, in 
ref. 0 , rapid rotation was added to a spherically symmetric 
solution and the large pressure depletion inevitably produced a 
dependence of the results on the amount of pressure reduction. 
Hereafter we will concentrate on the gravitational-wave emis¬ 
sion, but a discussion of the dynamics of matter and trapped 
surfaces can be found in 0 . 

While several different methods are possible for the extrac¬ 
tion of the gravitational-radiation content in numerical space- 
times, we have adopted a gauge-invariant approach in which 
the spacetime is matched with the non-spherical perturbations 
of a Schwarzschild black hole (see refs. 0 IT(1 JTl] for ap¬ 
plications to Cartesian coordinates grids). In practice, a set 
of “observers” is placed on 2 -spheres of fixed coordinate ra¬ 
dius r ex , where they extract the gauge-invariant, odd and 
even-parity metric perturbations j 12tl . Here £, m are the 
indices of the angular decomposition and we usually compute 
modes up to £ = 5 with m = 0; modes with m ^ 0 are essen¬ 
tially zero because of the high degree of axisymmetry in the 
collapse. Validations of this approach in 3D vacuum space- 
times can be found in refs. fllOt Il3l Il4ll. while its use with 
matter sources has first been reported in 1 1511 . 

Although the position of such observers is arbitrary and the 
information they record must be the same for waves extracted 
in the wave-zone, we place our observers between 40 M and 
50 M from the centre of the grid so as to maximize the length 
of the extracted waveform. We note that while a similar choice 
was made in ref. 0, it still provides only an approximate de¬ 
scription of what would be observed at spatial infinity. 

Using the odd and even-parity perturbations Qf m = A 
and Q\ m = AT'^, where A = ^2(£ + 2 )\/{£ - 2)!, we re¬ 
port in Fig. Qtlic lowest-order multipoles for Qj rn with the off¬ 
set produced by the stellar quadrupole removed 0 . The left 
panel, in particular, refers to the £ = 2 mode as extracted by 
four different observers at increasing distances and expressed 
in retarded time. The right panel instead refers to the £ = 4 
mode, with the inset giving a comparison between the two 
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FIG. 1: The left panel shows the 1 = 2, even-parity perturbation as extracted by observers at different positions r ex expressed in retarded time. 
The right panel refers to the l = 4 mode, with the inset offering a comparison in the amplitudes of the two modes. 


modes and showing that the gravitational-wave signal is es¬ 
sentially quadrupolar, with the 1 = 2 mode being about an 
order of magnitude larger than the £ = 4 mode. 

The very good overlap of the waveforms measured at dif¬ 
ferent positions is important evidence that the extraction has 
been performed in the wave-zone, since the invariance under 
a retarded-time scaling is a property of the solutions of a wave 
equation. The overlap disappears if the outer boundary is too 
close or when the waves are extracted at smaller radii. A sim¬ 
ilar overlap is seen also for the £ = 4 mode (not shown). 

Another indication that the waveforms in Fig. ^ are an ac¬ 
curate description of the gravitational radiation produced by 
the collapse comes by analysing their power spectra. The 
collapse, in fact, can be viewed as the rapid transition be¬ 
tween the spacetime of the initial equilibrium star and the 
spacetime of the produced rotating black hole. It is natu¬ 
ral to expect, therefore, that the waveforms produced in this 
process will reflect the basic properties of both spacetimes 
and in particular the fundamental frequencies of oscillation. 
We validate this in Fig. |3 where we show the power spec¬ 
tral densities (PSD) of the waveforms of the metric perturba¬ 
tions Q j, and Q^ 0 reported in Fig. ^ (the units on the (/-axis 
are arbitrary). The upper panel of Fig. [3 in particular, shows 
the PSD of Q 2 q and compares it with the frequencies of the 
£ = 2, to = 0 quasi-normal mode (QNM) of a Kerr black hole 
with M = 1.861 Mq and a = 0.6 |Hi] (dashed line at 6.7 
kHz) as well as with the first w u “interface” mode for a typi¬ 
cal compact star with M = 1.27 Mq and R = 8.86 km (dot¬ 
ted line at 8.8 kHz) El. Similarly, the lower panel of Fig.d 
shows the PSD of Q XIJ comparing it with the £ = 4, m = 0 
QNM of a Schwarzschild black hole Cl (dashed line at 14.0 
kHz) and the first W\ “curvature” mode 11711 (dotted line at 
12.8 kHz) Cl . It is quite apparent that both peaks in the PSDs 
are rather far from the maximum sensitivity area of modern 
interferometric and bar detectors. 


Although the waveforms have very short duration with very 
broad PSDs, Fig. |2] shows that these have strong and narrow 
peaks (a similar behaviour can be shown to be present also for 
the odd-parity modes d)- Indeed, the excellent agreement 
between the position of these peaks and the fundamental fre¬ 
quencies of the vacuum and non-vacuum spacetimes is an im¬ 
portant confirmation of the robustness of the results obtained. 

Using the extracted gauge-invariant quantities it is also pos¬ 
sible to calculate the transverse traceless (TT) gravitational- 
wave amplitudes in the two polarizations h + and h x as 

h +- ih -=^J2(Qtra- i f'QUt'yiA (1) 

r i,m ' ' 

where _ 2 Y tm is the s = —2 spin-weighted spherical har¬ 
monic. Because of the small amplitude of higher-order 
modes, the TT wave amplitudes can be simply expressed 
as h+ ~ h + (Q2 0 ,Qi t 0 ) and h x ~ h x (Q^ 0 ,Q^ 0 ), where 
Q 3 0 Qqq Their waveforms are shown in Fig. 0 for the 
detector at r ex = 37.1 M and for two different inclination 
angles d- Note that the amplitudes in the cross polariza¬ 
tion are about one order of magnitude smaller than those in 
the plus polarization, with the maximum amplitudes in a ra¬ 
tio \(r/M)h x | max /|(r/M)/i + | max ~ 0.06. Here, the odd- 
parity perturbations, which are zero in spacetime with axial 
and equatorial symmetries, are just the result of the coupling, 
induced by the rotation, with the even-parity perturbations. 

A precise comparison of the amplitudes in Fig. [3] with the 
corresponding ones calculated in |lj] is made difficult by the 
differences in the choice of initial data and, in particular, by 
the impossibility of reaching a > 0.54 when modelling con¬ 
sistently stationary polytropes in uniform rotation. However, 
when interpolating the results in d for the relevant values of 
a, we find a very good agreement in the form of the waves, 
but also that our estimates are about one order of magnitude 
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FIG. 2: Power spectra of the waveforms reported in Fig. Q The 
dashed vertical lines indicate the frequencies of the QNMs of a black 
hole, while the dotted ones the w\ and w n modes of a typical star. 

smaller, with \(r/M)h+\ max ~ 0.00225. Furthermore, we 
observe the amplitude of the gravitational waves to increase 
with the pressure reduction, thus suggesting that the origin of 
the difference is related mainly to this j^]. 

Following ref. IH and considering the optimal sensitiv¬ 
ity of VIRGO for the burst signal only, we set an upper limit 
for the characteristic amplitude produced in the collapse of 
a rapidly and uniformly rotating poly tropic star at 10 kpc to 
be h c = 5.77 x 10 _22 (M/M Q ) at a characteristic frequency 
f c = 931 Hz. In the case of LIGO I, instead, we obtain 
h c = 5.46 x 10 ~ 22 (M/Mq) at f c = 531 Hz. In both cases, 
the signal-to-noise ratio is S/N ~ 0.25, but this can grow to 
be < 4 in the case of LIGO II. These ratios could be increased 
considerably with the detection of the black hole ringing fol¬ 
lowing the initial burst @1. Computing the emitted power as 


dE 

dt 


1 

327T 


E 


dQL 

dt 



( 2 ) 


the total energy lost to gravitational radiation is E = 1.45 x 
10 ~ 6 (M/Mq). This is about two orders of magnitude smaller 
than the estimate made in |l|] for a star with a = 0.54, but 
larger than the energy losses computed recently in the collapse 
of rotating stellar cores to protoneutron stars |2£)]. 

In conclusion, we have presented the first waveforms from 
the gravitational collapse of rapidly rotating stars to black 
holes using 3D grids with Cartesian coordinates. The great 
potential shown by the PMR techniques employed here opens 
the way to a number of applications that would be otherwise 
intractable with uniform Cartesian grids. Work is in progress 
to consider initial models with realistic EOSs or in differen¬ 
tial rotation, for which values a > 1 can be reached and more 
intense gravitational radiation is expected. 
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FIG. 3: TT gravitational-wave amplitudes extracted by the detector 
at r ex = 37.1 M and for two different inclination angles. 
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